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The R136 star cluster hosts several stars whose individual masses 
greatly exceed the accepted 150 M© stellar mass limit 
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ABSTRACT 

Spectroscopic analyses of hydrogen-rich WN5-6 stars within the young star clusters NGC 
O ■ 3603 and R136 are presented, using archival Hubble Space Telescope and Very Large 

Telescope spectroscopy, and high spatial resolution near-IR photometry, including Multi- 
^ . Conjugate Adaptive Optics Demonstrator (MAD) imaging of R136. We derive high stellar 

f^" temperatures for the WN stars in NGC 3603 (T* - 42±2 kK) and R136 (T* - 53± 3 kK) 

plus dumping-corrected mass-loss rates of 2-5 xlO~^ M© yr~^ which closely agree with 
theoretical predictions from Vink et al. These stars make a disproportionate contribution to 
K^ I the global ionizing and mechanical wind power budget of their host clusters. Indeed, R136al 

"^" . alone supplies r^l% of the ionizing flux of the entire 30 Doradus region. Comparisons with 

2^ ■ stellar models calculated for the main- sequence evolution of 85 - 500 M© accounting for 

rotation suggest ages of ~ 1.5 Myr and initial masses in the range 105 - 170 M© for three 
systems in NGC 3603, plus 165 - 320 M© for four stars in R136. Our high stellar masses 
r^ ' are supported by consistent spectroscopic and dynamical mass determinations for the com- 

^^ \ ponents of NGC 3603 Al. We consider the predicted X-ray luminosity of the R136 stars if 

^p ■ they were close, colliding wind binaries. R136c is consistent with a colliding wind binary 

system. However, short period, colliding wind systems are excluded for R136a WN stars if 
mass ratios are of order unity. Widely separated systems would have been expected to harden 
owing to early dynamical encounters with other massive stars within such a high density en- 
vironment. From simulated star clusters, whose constituents are randomly sampled from the 
Kroupa initial mass function, both NGC 3603 and R136 are consistent with an tentative upper 
mass limit of ~300 M©. The Arches cluster is either too old to be used to diagnose the upper 
mass limit, exhibits a deficiency of very massive stars, or more likely stellar masses have been 
underestimated - initial masses for the most luminous stars in the Arches cluster approach 200 
M0 according to contemporary stellar and photometric results. The potential for stars greatly 
exceeding 150 M© within metal-poor galaxies suggests that such pair-instability supernovae 
could occur within the local universe, as has been claimed for SN 2007bi. 

Key words: binaries: general - stars: early-type - stars: fundamental parameters - stars: 
Wolf-Rayet - galaxies: star clusters: individual (R136, NGC 3603, Arches) 



1 INTRODUCTION of a limit close to 150 M© (Figer 2005; Koen 2006). If stars above 

^ , , , , ^ . ^ . . , , this limit were to exist, they would be exclusive to the youngest, 

Although the formation 01 very massive stars remains an unsolved , • , . . 1 . , • , . /t- ^^^cn 

^ ^ ^ , . ,^. , n xr 1 ^^^^x , t t highest mass star clusters, which are very compact (Figer 2005). 

problem 01 astrophysics (Zmnecker & Yorke 2007), the past decade ^t .• h 1 j • • r 1 1 . • .1 

f ,.pp ,,r, , • 1 1 However, spatially resolved imaging 01 such clusters is currently 

has seen a shiit irom a beliei that there is no observational upper 1 • . .1 A/rn ^.^t j •. . n-. 1 • 1 .1 

^^, ^^ ^^^ox 1 -t t exclusive to the Milky Way and Its satellite galaxies, where they are 

stellar mass cutoii (e.s.Massey 2003) to the widespread acceptance t jj-.- .1 . j. • .• r . n 

^ -^ ^ ^ very rare. In addition, the accurate determination 01 stellar masses 

generally relies upon spectroscopic and evolutionary models, un- 
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less the star is a member of an eclipsing binary system (Moffat 
2008), something which occurs extremely rarely. 

Most, if not all, stars form in groups or clusters (Lada & Lada 
2003). An average star forms with an initial mass of ^0.5 M© 
while the relative proportion of stars of higher and lower mass 
obeys an apparently universal initial mass function (IMF, Kroupa 
2002). In addition, there appears to be a relationship between the 
mass of a cluster and its highest-mass star (Weidner & Kroupa 
2006, Weidner et al. 2010). High mass stars (>8 M©), ultimately 
leading to core-collapse supernovae, usually form in clusters ex- 
ceeding 100 M while stars approaching 150 M© have been de- 
tected in two 10^ M© Milky Way clusters, namely the Arches 
(Figer 2005; Martins et al. 2008) and NGC 3603 (Schnurr et al. 
2008a). Is this 150 M© limit statistical or physical? 

The determination of stellar masses of very massive stars from 
colour-magnitude diagrams is highly unreliable, plus corrections to 
present mass estimates need to be applied to estimate initial masses. 
Studies are further hindered by severe spatial crowding within the 
cores of these star clusters (Maiz Apellaniz 2008). In addition, so- 
phisticated techniques are required to extract the physical parame- 
ters of early-type stars possessing powerful stellar winds from opti- 
cal/infrared spectroscopy (Conti et al. 2008). In general, the incor- 
poration of line blanketing has led to a reduction in derived temper- 
atures for O stars from photospheric lines (Puis et al. 2008) while 
the reverse is true for emission lines in Wolf-Rayet stars (Crowther 
2007). 

Here we re-analyse the brightest members of the star cluster 
(HD 97950) responsible for NGC 3603, motivated in part by the 
case of Al which is an eclipsing binary system, whose individual 
components have been measured by Schnurr et al. (2008a). We also 
re-analyse the brightest members of R136 (HD 38268) - the cen- 
tral ionizing cluster of the Tarantula nebula (30 Doradus) within 
the Large Magellanic Cloud (LMC). R136 is sufficiently young 
(Massey & Hunter 1998) and massive (<5.5 x 10^ M©, Hunter 
et al. 1995) within the Local Group of galaxies to investigate the 
possibility of a physical limit beyond 150 M© (Massey & Hunter 
1998, Selman et al. 1999; Oey & Clarke 2005). 

As recently as the 1980s, component 'a' within R136 was be- 
lieved to be a single star with a mass of several thousand solar 
masses (Cassinelli et al. 1981; Savage et al. 1983), although others 
favoured a dense star cluster (Moffat & Seggewiss 1983). The lat- 
ter scenario was supported by speckle interferometric observations 
(Weigelt & Bauer 1985) and subsequently confirmed by Hubble 
Space Telescope (HST) imaging (Hunter et al. 1995). 

Here, we present new analyses of the brightest sources of 
NGC 3603 and R136. In contrast to classical Wolf-Rayet stars, 
these WN stars are believed to be young, main-sequence stars, al- 
beit possessing very strong stellar winds as a result of their high 
stellar luminosities (Crowther 2007). Section 2 describes archival 
HST and Very Large Telescope (VLT) observations, while the spec- 
troscopic analysis is presented in Sect 3. Comparisons with con- 
temporary evolutionary models are presented in Sect. 4, reveal- 
ing spectroscopic masses in excellent agreement with dynamical 
masses of the components of NGC 3603 Al. For the R136 stars, 
exceptionally high initial masses of 165 - 320 M© are inferred. 
We consider the possibility that these stars are binaries in Sect 5, 
from archival X-ray observations of R136. In Sect. 6 we simulate 
star clusters to re-evaluate the stellar upper mass limit, and find that 
both R136 and NGC 3603 are consistent with a tentative upper limit 
of ~ 3OOM0. Sect. 7 considers the contribution of the R136 stars 
to the global properties of both this cluster and 30 Doradus. Finally, 



Table 1. Log of spectroscopic observations of R136 and NGC 3603 stars 
used in this study 



Star 


Instrument 


Grating 


Date 


Proposal/PI 


R136al 


HST/HRS 


G140L 


1994 Jul 


5157/Ebbets 


R136a3 


HST/HRS 


G140L 


1994 Jul 


5157/Ebbets 


R136al 


HST/FOS 


G400H, 
G570H 


1996 Jan 


6018/Heap 


R136a2 


HST/FOS 


G400H, 
G570H 


1996 Jan 


6018/Heap 


R136a3 


HST/FOS 


G400H, 
G570H 


1996 Jan 


6018/Heap 


R136C 


HST/FOS 


G400H 


1996 Nov 


6417/Massey 


R136al+a2 


VLT/SINFONI 


K 


2005 Nov-, 
2005 Dec 


076.D-0563/Schnurr 


R136a3 


VLT/SINFONI 


K 


2005 Nov- 
2005 Dec 


076.D-0563/Schnurr 


R136C 


VLT/SINFONI 


K 


2005 Nov- 
2005 Dec 


076.D-0563/Schnurr 


NGC 3603A1 


HST/FOS 


G400H 


1994 Sep 


5445/Drissen 


NGC 3603B 


HST/FOS 


G400H 


1994 Sep 


5445/Drissen 


NGC 3603C 


HST/FOS 


G400H 


1994 Sep 


5445/Drissen 


NGC 3603A1 


VLT/SINFONI 


K 


2005 Apr- 

2006 Feb 


075.D-0577/Moffat 


NGC 3603B 


VLT/SINFONI 


K 


2005 Apr- 

2006 Feb 


075.D-0577/Moffat 


NGC 3603C 


VLT/SINFONI 


K 


2005 Apr- 

2006 Feb 


075.D-0577/Moffat 



Table 2. Stars brighter than Mks ~ -6 mag within 13 arcsec (0.5 parsec) of 
NGC 3603 Al, together with V-band photometry from Melena et al. (2008) 
in parenthesis. 



Name Sp Type 



(my) 



(Ay) 



M^^ Binary^ 

(My) 



Al WN6h 7.42 ±0.05 0.59 ± 0.03 -7.57 ± 0.12 Yes 

(11.18) (4.91±0.25) (-8.13±0.27) 

B WN6h 7.42 ±0.05 0.56 ± 0.03 -7.54 ± 0.12 No? 

(11.33) (4.70±0.25) (-7.77±0.27) 

C WN6h^ 8.28 ±0.05 0.56 ± 0.03 -6.68 ±0.12 Yes 

(11.89) (4.66±0.25) (-7.17±0.27) 

(a) For a distance of 7.6 ± 0.35 kpc (distance modulus 14.4 ±0.1 mag) 

(b) Al is a 3.77 day double-eclipsing system, while C is a 8.9 day SBl 
(Schnurr et al. 2008a) 

(c) An updated classification scheme for Of, Of/WN and WN stars (N.R. 
Walborn, & R A. Crowther, in preparation) favours 03 If*/WN6 for NGC 
3603 C 



we consider the broader significance of an increased mass limit for 
stars in Sect. 8. 



2 OBSERVATIONS 

Our analysis of WN stars in NGC 3603 and R136 is based upon 
archival UV/optical HST and near-IR VLT spectroscopy, sum- 
marised in Table [T] combined with archival high spatial resolution 
near-IR imaging. We prefer the latter to optical imaging due to re- 
duced extinction, efficient correction for severe spatial crowding 
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Figure 1. VLT MAD Ks-band 12 x 12 arcsec (3x3 parsec for the LMC 
distance of 49 kpc) image of R 136 (Campbell et al. (2010) together with a 
view of the central 4x4 arcsec (1x1 parsec) in which the very massive 
WN5h stars discussed in this letter are labelled (component b is a lower 
mass WN9h star). Relative photometry agrees closely with integral field 
SINFONI observations (Schnurr et al. 2009). 



(Schnurr et al. 2008a, 2009) and consistency with studies of other 
young, high mass clusters (e.g. Arches, Martins et al. 2008). 



2.1 NGC 3603 

Our spectroscopic analysis of the three WN6h systems within NGC 
3603 is based upon archival HST/FOS spectroscopy from Dris- 
sen et al. (1995) plus integral field VLT/SINFONI near-IR spec- 
troscopy from Schnurr et al. (2008a) - see Table [T] for the log of 
observations. The spectral resolution of the SINFONI datasets is 
R^3000, with adaptive optics (AO) used to observe Al, B and C, 
versus R'^1300 for FOS for which a circular aperture of diameter 
0.26 arcsec was used. 

Differential photometry from VLT/SINFONI integral field 
observations (Schnurr et al 2008a), were tied to unpublished 
VLT/ISAAC Ks-band acquisition images from 16 Jun 2002 (ESO 
Programme 69.D-0284(A), P.I. Crowther) using the relatively iso- 
lated star NGC 3603 C. The VLT/ISAAC frames are calibrated 
against 2MASS photometry (Skrutskie et al. 2006) using 6 stars 
in common (±0.05 mag). These are presented in Table |2] together 
with absolute magnitudes resulting from a distance of 7.6±0.35 
kpc (distance modulus 14.4db0.1 mag) plus an extinction law from 
Melena et al. (2008) using 



A^ 



1.17?^"^ + (0.29 



0.35)i?^^^^^°^ 



where Ry^ = 3.1 for the Milky Way foreground component and 
j^NGC3603 ^ 4 3 f^j. ^j^g internal NGC 3603 component. We obtain 
Axg = 0.12 Av ^ 0.56 - 0.59 mag from spectral energy distribu- 
tion fits, which are consistent with recent determinations of E(B-V) 
= 1.39 mag from Melena et al. (2008). In their analysis, Crowther 
& Dessart (1998) used a lower overall extinction of E(B-V) = 1.23 
mag, albeit a higher distance of 10 kpc (distance modulus of 15.0 
mag) to NGC 3603. 



2.2 R136 

Our UV/optical/infrared spectroscopic analysis of all four 
hydrogen-rich WN5h Wolf-Rayet (WR) stars in R136 (Crowther 
& Dessart 1998) is also based upon archival HST and VLT spec- 
troscopy - see Table [T]for the log of observations. 

Goddard High Resolution Spectrograph (GHRS) ultraviolet 
observations of R 136a stars have been described by de Koter et 
al. (1997), achieving a spectral resolution of R^ 2000-3 000 using 
the SSA aperture of 0.22x0.22 arcsec. This prevented individual 
UV spectroscopy for R136al and a2 which are separated by ^0.1 
arcsec. For context, 0.1 arcsec subtends 5000 AU at the distance 
of the LMC. R136a2 UV spectroscopy was attempted by S. Heap 
(Programme 5297), but was assessed to be unreliable by de Koter 
et al. (1997) and so is also excluded here. Note also that R136c has 
not been observed with GHRS. Visual Faint Object Spectroscopy 
(FOS) datasets, with a circular aperture of diameter 0.26 arcsec 
achieved R^1300, although R136al and R136a2 once again suffer 
significant contamination from one another. It is solely at near-IR 
wavelengths (SINFONI) that R136al and a2 are spectrally sepa- 
rated, for which R136b served as an AO reference star (Schnurr 
et al. 2009). The spectral resolution of the SINFONI datasets is 
R-3000. 

We employ high spatial resolution Ks-band photome- 
try of R136. Differential Ks photometry from AO assisted 
VLT/SINFONI integral field datasets are tied to identical spa- 
tial resolution wider field VLT Multi-Conjugate Adaptive Optics 
Demonstrator (MAD) imaging (Campbell et al. 2010) using the rel- 
atively isolated star R136b (WN9h). Three overlapping Fields were 
observed with VLT/MAD, for which Field 1 provided the highest 
quality in R136 (FWHM ^-^ 0.1 arcsec), as shown in Fig. [T] itself 
calibrated using archival HAWK-I and 2MASS datasets (see Camp- 
bell et al. 2010). For 12 stars in common between MAD photom- 
etry and HST Near Infrared Camera and Multi Object Spectrome- 
ter (NICMOS) F205W imaging (Brandner et al. 2001) transformed 
into the CTIO K-band system, Campbell et al. (2010) find mx^ 
(MAD) - mK (HST) = - 0.04 ±0.05 mag. 

In Table [3] we present Ks aperture photometry and inferred ab- 
solute magnitudes of stars brighter than Mk^ "^ -6 mag within 20 
arcsec (5 parsec) of R136al. Of these only R136a, b and c compo- 
nents lie within a projected distance of 1 pc from R136al. Spec- 
tral types are taken from Crowther & Dessart (1998) or Walborn & 
Blades (1997), although other authors prefer alternative nomencla- 
ture, e.g. 02 If for both Mk39 and R136a5 according to Massey et 
al. (2004, 2005). 

Interstellar extinctions for the R136 WN5h stars are derived 
from UV to near-IR spectral energy distribution fits (R136c lacks 
UV spectroscopy), adopting foreground Milky Way (LMC) extinc- 
tions of Aks = 0.025 (0.06) mag, plus variable internal 30 Doradus 
nebular extinction. The adopted extinction law follows Fitzpatrick 
& Savage (1984) as follows 



Av = 0.07Rt 



+ 0.16i?v^^ 



+ (0.25 ■ 



■ 0.45)i???°°^ 



where R^^ = R^^^ = 3.2 and Rf °°^ = 4.0. We derive Ak, -0.22 
mag for the R136a stars and 0.30 mag for R136c and adopt Ak^ = 
0.20 mag for other stars except that Mk34 (WN5h) mirrors the 
higher extinction of R136c. From Table[3] R136c is 0.5 mag fainter 
than R136a2 in the V-band (Hunter et al. 1995) but is 0.06 mag 
brighter in the K^-band, justifying the higher extinction. An analy- 
sis based solely upon optical photometry could potentially underes- 
timate the bolometric magnitude for R136c with respect to R136a2. 
The main source of uncertainty in absolute magnitude results from 
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Figure 2. Spectral energy distributions of 
NGC 3603 WN6h stars from HST/FOS 
together with Ks photometry from 
VLT/SINFONI and VLT/ISAAC plus 
reddened theoretical spectral energy 
distributions (red lines). Components Ala 
(green lines) and Alb (blue lines) are also 
shown separately. 



the distance to the LMC. The mean of 7 independent techniques 
(Gibson 2000) suggests an LMC distance modulus of 18.45 ± 0.06 
mag. Here we adopt an uncertainty of ±0.18 mag owing to system- 
atic inconsistencies between the various methods. 



3 SPECTROSCOPIC ANALYSIS 

For our spectroscopic study we employ the non-LTE atmosphere 
code CMFGEN (Hillier & Miller 1998) which solves the radia- 
tive transfer equation in the co-moving frame, under the additional 
constraints of statistical and radiative equilibrium. Since CMFGEN 
does not solve the momentum equation, a density or velocity struc- 
ture is required. For the supersonic part, the velocity is parameter- 
ized with an exponent of /3 = 0.8. This is connected to a hydro- 
static density structure at depth, such that the velocity and velocity 
gradient match at this interface. The subsonic velocity structure is 
defined by a fully line-blanketed, plane-parallel TLUSTY model 



(Lanz & Hubeny 2003) whose gravity is closest to that obtained 
from stellar masses derived using evolutionary models, namely 
log^ = 4.0 for R136 stars and log^ = 3.75 for NGC 3603 stars. 
CMFGEN incorporates line blanketing through a super-level ap- 
proximation, in which atomic levels of similar energies are grouped 
into a single super-level which is used to compute the atmospheric 
structure. 

Stellar temperatures, T*, correspond to a Rosseland optical 
depth 10, which is typically 1,000 K to 2,000 K higher than ef- 
fective temperatures T2/3 relating to optical depths of 2/3 in such 
stars. 

Our model atom include the following ions: H I, He I-II, C III- 
IV, N III-V, O III- VI, Ne IV- V, Si IV, P IV- V, S IV- V, Ar V-VII, Fe 
IV- VII, Ni V-VII, totalling 1,141 super-levels (29,032 lines). Other 
than H, He, CNO elements, we adopt solar abundances (Asplund 
et al. 2009) for NGC 3603, which are supported by nebular studies 
(Esteban et al. 2005, Lebouteiller et al. 2008). LMC nebular abun- 
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Table 3. Stars brighter than M^^ ~ -6 mag within 20 arcsec (5 parsec) 
of R136al, together with V-band photometry from Hunter et al. (1995) in 
parenthesis. 



Name 


Sp Type 


(my) (Av) 


(My) 


Binary^ 


R134 


WN6(h) 


10.91 ±0.09 0.21 ±0.04 
(12.89±0.08) (1.77±0.33) 


-7.75 ± 0.20 
(-7.33±0.38) 


No? 


R136al 


WN5h 


11.10 ±0.08 0.22 ±0.02 
(12.84±0.05) (1.80±0.17) 


-7.57 ± 0.20 
(-7.41±0.25) 


No? 


R136C 


WN5h 


11.34 ±0.08 0.30 ±0.02 
(13.47±0.08) (2.48±0.17) 


-7.41 ± 0.20 
(-7.46±0.26) 


Yes? 


R136a2 


WN5h 


11.40 ±0.08 0.23 ±0.02 
(12.96±0.05) (1.92±0.17) 


-7.28 ± 0.20 
(-7.41±0.25) 


No? 


Mk34 


WN5h 


11.68 ±0.08 0.27 ±0.04 
(13.30±0.06) (2.22±0.33) 


-7.04 ± 0.20 
(-7.37±0.38) 


Yes 


R136a3 


WN5h 


11.73 ±0.08 0.21 ±0.02 
(13.01±0.04) (1.72±0.17) 


-6.93 ± 0.20 
(-7.16±0.25) 


No? 


R136b 


WN9ha 


11.88 ±0.08 0.21 ±0.04 
(13.32±0.04) (1.74±0.33) 


-6.78 ± 0.20 
(-6.87±0.38) 


No? 


Mk39 


02-3 If AVN^ 


12.08 ±0.08 0.18 ±0.04 
(13.01±0.08) (1.46±0.33) 


-6.55 ± 0.20 
(-6.90±0.38) 


Yes 


Mk42 


02-3IfAVN^ 


12.19 ±0.08 0.17 ±0.04 
(12.84±0.05) (1.38±0.33) 


-6.43 ± 0.20 
(-6.99±0.38) 


No? 


Mk37a 


04If+^ 


12.39 ±0.11 0.21 ±0.04 
(13.57±0.05) (1.74±0.33) 


-6.27 ± 0.21 
(-6.62±0.38) 


No? 


Mk37Wa04If+ 


12.39 ±0.11 0.19 ±0.04 


-6.25 ± 0.21 


? 






(13.49±0.05) (1.62±0.33) 


(-6.58±0.38) 




R136a5 


02-3IfAVN^ 


12.66 ± 0.08 0.21 ± 0.04 
(13.93±0.04) (1.74±0.33) 


-6.00 ± 0.20 
(-6.26±0.38) 


No? 



(a) For a distance modulus of 18.45 ±0.18 mag (49 ± 4 kpc) 

(b) R136c is variable and X-ray bright, Mk39 is a 92-day SBl binary 
(Massey et al. 2002; Schnurr et al. 2008b) and Mk 34 is a binary accord- 
ing to unpublished Gemini observations (O. Schnurr et al. in preparation) 

(c) An updated classification scheme for Of, OfAVN and WN stars (N.R. 
Walborn, & R A. Crowther, in preparation) favours 02 If* for Mk 42 and 
R136a5, 03.5 If for Mk 37a and 02 IfAVN for Mk 39 

dances (Russell & Dopita 1990) are adopted for R136, with other 
metals scaled to 0.4 Z©, also supported by nebular studies of 30 
Doradus (e.g. Peimbert 2003, Lebouteiller et al. 2008). We have 
assumed a depth-independent Doppler profile for all lines when 
solving for the atmospheric structure in the co-moving frame, while 
in the observer's frame, we have adopted a uniform turbulence of 
50 km s~^. Incoherent electron scattering and Stark broadening for 
hydrogen and helium lines are adopted. With regard to wind clump- 
ing (Hillier 1991), this is incorporated using a radially-dependent 
volume filling factor, /, with /oo =0.1 at i;oo, resulting in a reduc- 
tion in mass-loss rate by a factor of y^(l//) ^ 3. 



3.1 NGC 3603 

Fig.[2]compares spectral energy distributions for the WN6h stars in 
NGC 3603 with reddened, theoretical models, including the indi- 
vidual components of AL For this system we obtain Mks = -7.57 
mag, such that we adopt Am = mAia - mAib = -0.43 ±0.3 mag 
for the individual components (Schnurr et al. 2008a). We have es- 
timated individual luminosities in two ways. First, on the basis of 
similar mean molecular weights /i for individual components, stel- 
lar luminosities result from L oc M"" , with a ^ 1.5 for ZAMS mod- 
els in excess of 85 M©. We derive L(Ala)/L(Alb) = 1.62+g;^ from 
the dynamical masses of 116 ± 31 M© and 89 ± 16 M© obtained 
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Figure 3. Rectified, spatially resolved near-IR (VLT/SINFONI) spec- 
troscopy of NGC 3603 WN6h stars (Schnurr et al. 2008a, black), includ- 
ing Al close to quadrature (Ala blueshifted by 330 km s"-*^ and Alb red- 
shifted by 433 km s~^), together with synthetic spectra (red, broadened by 
50 km s-i). N III 2.103 - 2.1 15//m and He I 2.112 - 2.1 13/im contribute to 
the emission feature blueward of Br 7 



by Schnurr et al. (2008a) for Ala and Alb, respectively. Alterna- 
tively, we have obtained a luminosity ratio for components within 
Al from a fit to the NICMOS lightcurve (Moffat et al. 2004, and 
E. Antokhina, priv. comm), using the radii of the Roche lobes and 
lightcurve derived temperatures. We also obtain L(Ala)/L(Alb) = 
1.62 from this approach. 

Although we have matched synthetic spectra to near-IR pho- 
tometry, we note that no major differences would be obtained with 
either HST Wide Field and Planetary Camera 2 (WFPC2) or Ad- 
vanced Camera for Surveys (ACS) datasets. For example, our red- 
dened spectral energy distribution for NGC 3603B implies B = 
12.40 mag, which matches HST/ACS photometry to within 0.06 
mag (Melena et al. 2008). 

We have used diagnostic optical and near-IR lines, includ- 
ing Nlll 4634-41, 2.103-2.115/im, Niv 3478-83, 4058, together 
with Hell 4686, 2.189/im plus Br7. We are unable to employ he- 
lium temperature diagnostics since He I lines are extremely weak 
at optical and near-IR wavelengths. Spectroscopic comparisons 
with HST/FOS datasets are presented in Fig. |2] Overall, emis- 
sion features are well reproduced, although absorption compo- 
nents of higher B aimer-Pickering lines are too strong, especially 
for NGC 3603C. Fits are similar to Crowther & Dessart (1998), in 
spite of an improved TLUSTY structure within the photosphere. 
Fortunately, our spectral diagnostics are not especially sensitive to 
the details of the photosphere since they assess the inner wind con- 
ditions. Spectroscopically, we use VLT/SINFONI spectroscopy of 
Ala and Alb obtained close to quadrature (Schnurr et al. 2008a) 
to derive mass-loss rates and hydrogen contents. The ratio of He II 
2.189/im to Br7 provides an excellent diagnostic for the hydrogen 
content in WN stars, except for the latest subtypes (approximately 
WN8 and later). As such, we are able to use primarily optical spec- 
troscopy for the determination of stellar temperatures, with hydro- 
gen content obtained from near-IR spectroscopy. Near-IR compar- 
isons between synthetic spectra and observations are presented for 
each of the WN6h stars within NGC 3603 in Fig. [3] 

Physical properties for these stars in NGC 3603 are presented 
in Table ID including evolutionary masses obtained from solar- 
metallicity non-rotating models. Errors quoted in the Table account 
for both photometric and spectroscopic uncertainties, and repre- 
sent the range of permitted values. However, for NGC 3603 C we 
assume that the primary dominates the systemic light since the 



p. A. Crowther et al. 




R136al 
Av=1.8 



1.6 

^1.2 



''■r''^''^'*^--^ VV_A 





4000 4250 4500 4750 



log X (A) 



Figure 4. Spectral energy distributions of 
R136 WN5h stars from HST/FOS together 
using Ks photometry from VLT/SINFONI 
calibrated with VLT/MAD imaging. Red- 
dened theoretical spectral energy distribu- 
tions are shown as red lines. 



companion is not detected in SINFONI spectroscopy (Schnurr et 
al 2008a). Should the companion make a non-negligible contribu- 
tion to the integrated light, the derived properties set out in Table |4] 
would need to be corrected accordingly. The improved allowance 
for metal line blanketing implies ^10-20% higher stellar temper- 
atures (T* -- 40,000 - 44,000 K) and, in turn, larger bolometric 
corrections (MboI - My '^ -3.7 db 0.2 mag) for these stars than 
earlier studies (Crowther & Dessart 1998), who adopted a different 
combination of E(B-V) and distance for NGC 3603. We note that 
Schmutz & Drissen (1999) have previously derived T* ^ 46,000 
K for NGC 3603 B, resulting in a similar luminosity (logL/L© 
= 6.4) to that obtained here. Analysis of the NICMOS lightcurve 
(Moffat et al. 2004) yields T(Ala)/T(Alb) = 1.06, in support of 
the spectroscopically-derived temperature ratio, albeit 7% higher 
in absolute terms. In view of the underlying assumptions of the 
photometric and spectroscopic techniques, overall consistency is 
satisfactory. 



For Ala and B, we adopt abundances of Xc = 0.008%, Xn 
= 0.8% and Xo = 0.013% by mass, as predicted by evolutionary 
models for Xh -- 60%, versus Xc = 0.005%, Xn = 0.6% and Xo 
= 0.25% for Alb and C, for which Xh -^ 70%. 



3.2 R136 

Theoretical spectral energy distributions of all four WN5h stars are 
presented in Fig.lU Although we have matched synthetic spectra to 
VLT/MAD + SINFONI photometry, we note that no significant dif- 
ferences would be obtained from WFPC2 imaging (e.g. Hunter et 
al. 1995), providing these are sufficiently isolated. For example, our 
reddened spectral energy distribution for R136a3 implies V = 13.0 
mag, in agreement with F555W photometry. VLT/MAD photome- 
try is preferred to WFPC2 for the very crowded core of R136 (e.g. 
al and a2), although recent Wide Field Camera 3 (WFC3) imaging 
of R136 achieves a similar spatial resolution at visible wavelengths. 
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Table 4. Physical Properties of NGC 3603 WN6h stars. 



Table 5. Physical Properties of R136 WN5h stars. 



Name 


Ala 


Alb 


B 


C 


Name 
BAT99 


al 
108 


a2 
109 


a3 
106 


c 
112 


T* (kK)« 


42 ±2 


40 ±2 


42 ±2 


44 ±2 


log (L/Lo) 


6.39 ±0.14 


6.18 ±0.14 


6.46 ± 0.07 


6.35 ± 0.07 


T* (kK)« 


53 ±3 


53 ±3 


53 ±3 


51±5 


Rr = 2/3 (R©) 29.41^3^ 
NLyC (10^0 ,-1) 1,6+0.8 

M(lO-5M0yr-i) S.l^li 
logM-logM^.^i^ +0.14 


2CQ+7.2 
ZD.y_3 -^ 

0-85+S-23 

1 Q+0.9 
^•^-0.4 

+0.24 


33.8+2-^ 

1 Q+O^S 
^•^-0.3 

c 1+0.6 
-^•^-0.6 

+0.22 


26 2+21 

1-^-0.3 
1-^-0.2 

-0.04 


log(L/L0) 

Rt = 2/3 (Rq) 
NLyC (10^0 S-1) 

M(lO-5M0yr- 


6.94 ± 0.09 

^S 4+4-0 

6.6111 

-1^ c 1+0.9 
) -^-1-0. 8 


6.78 ± 0.09 
29.5tll 

A Q + 0.8 
4 6+0-8 


6.58 ± 0.09 

23.41-1 

3.oi°-i 

3 7+0.7 


6.75 ±0.11 
30.61^.? 

4.2lg:^ 

4 5+1-0 

^•-^-0.8 


Voo (kms-1) 


2600 ± 150 


2600 ± 150 


2300 ± 150 


2600 ± 150 


logM-logM^.^^ +0.09 


+0.12 


+0.18 


+0.06 


Xif (%) 


60 ±5 


70 ±5 


60 ±5 


70 ±5 


Voo(kms-i) 


2600 ± 150 


2450 ± 150 


2200 ± 150 


1950 ± 150 


Mcurrent (M©)^ 

Mk, (mag)^ 


148+^0 

UOtTr 
-7.0 ± 0.3 


92tll 
-6.6 ± 0.3 


132tll 
-7.5 ±0.1 


137111 
-6.7 ±0.1 


Xh (%) 
Minit (Mo)^ 

Mcurrent (M©)^ 


40 ±5 

320+^go 

265tll 

-7.6 ± 0.2 


35 ±5 

240ttl 

I95tll 

-13 ± 0.2 


40 ±5 
1651- 

1351^0 
-6.9 ± 0.2 


30 ±5 

220+^^ 

1751^5 
-7.4 ± 0.2 












Mk^ (mag) 



(a) Corresponds to the radius at a Rosseland optical depth of tross = 10 

(b) Component C is a 8.9 day period SBl system (Schnurr et al. 2008a) 

(c) dM/dtvink relates to Vink et al. (2001) mass-loss rates for Z = Z© 

(d) Mks = -7.57 ± 0.12 mag for Al, for which we adopt Am = mAia 
- mAib = -0.43 ± 0.30 mag (Schnurr et al. 2008a). The ratio of their lu- 
minosities follows from their dynamical mass ratios together with L oc 
fiM^-^ (and is supported by NICMOS photometry from Moffat et al. 
2004). 



Nevertheless, we adhere to VLT/MAD + SINFONI for our primary 
photometric reference since 

(i) Solely VLT/SINFONI spectrally resolves R136al from 
R136a2 (Schnurr et al. 2009). We also note the consistent line to 
continuum ratios between optical and near-IR diagnostics, ruling 
out potential late-type contaminants for the latter; 

(ii) Uncertainties in dust extinction at Ks are typically 0.02 mag 
versus ^-^0.2 mag in the V-band, recalling Aks ^ 0.12^^. R136c 
is significantly fainter than R136a2 at optical wavelengths yet is 
brighter at Ks and possesses a higher bolometric magnitude (con- 
trast these results with the optical study of Riihling 2008); 

(iii) We seek to follow a consistent approach to recent studies of 
the Arches cluster which necessarily focused upon the Ks-band for 
photometry and spectroscopy (e.g. Martins et al. 2008). 

We estimate terminal wind velocities from optical and near- 
IR helium lines. Ultraviolet HST/GHRS spectroscopy has been ob- 
tained for the R136a stars (Heap et al. 1994; de Koter et al. 1997), 
although the severe crowding within this region results in the spec- 
tra representing blends of individual components for R136al and 
a2 (apparent in Fig.|4l). As a result, we favour velocities from FOS 
and SINFONI spectroscopy for these cases. 

In Fig. [5] we present UV/optical/near-IR spectral fits for the 
WN5h star R136a3, which is sufficiently isolated that contamina- 
tion from other members of R136a is minimal in each spectral re- 
gion (in contrast to R136al and a2 at UV/optical wavelengths). 
Consistent optical and near-IR fits demonstrate that no underly- 
ing red sources contribute significantly to the K-band SINFONI 
datasets. Diagnostics include Ov 1371, Sv 1501, Niv 3478-83, 
4058, Nv 4603-20, 2.10/im together with Hell 4686, 2.189/im 
plus Ha, H/3, Br7. For T* < 50,000 K both Ov 1371 and Nv 
2.11/im are underestimated for the R136a WN5h stars, while for 
T* > 56,000 K, NiV 3478-83 and 4057 become too weak, and 
S V 1501 becomes too strong. Therefore, we favour T* ^^53,000 
± 3,000 K, except that T* - 51,000 ±5,000 K is preferred for 
R136c since UV spectroscopy is not available and N V 2.100/im is 



(a) Corresponds to the radius at a Rosseland optical depth of tross = 10 

(b) Component R136c is probably a coUiding-wind massive binary. For a 
mass ratio of unity, initial (current) masses of each component would cor- 
respond to ~160 M© (~130 Mo) 

(c) dM/dtvink relates to Vink et al. (2001) mass-loss rates for Z = 0.43 Z© 
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Figure 5. Rectified, ultraviolet (HST/GHRS), visual (HST/FOS) and near- 
IR (VLT/SINFONI) spectroscopy of the WN5h star R136a3 together with 
synthetic UV, optical and near-infrared spectra, for T* = 50,000 K (red) 
and T* =55,000 K (blue). Instrumental broadening is accounted for, plus an 
additional rotational broadening of 200 km s~ ^ . 
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Figure 6. Rectified, spatially resolved near-IR (VLT/SINFONI, Schnurr et 
al. 2009) spectroscopy of R136 WN5h stars (black) together with synthetic 
spectra, accounting for instrumental broadening (100 kms~^) plus rota- 
tional broadening of 200 km s~^ for R136a2, a3 and c. Consistent hydrogen 
contents are obtained from the peak intensity ratio of Br7/HeII 2.189/im 
and optical (Pickering-B aimer series) diagnostics. Nitrogen emission in- 
cludes Nv 2.1 00//m and N III 2.103 - 2.115/im. 



weak/absent. Again, we are unable to employ helium temperature 
diagnostics since He I lines are extremely weak. 

Fig. m compares synthetic spectra with HST/FOS spec- 
troscopy, for which emission features are again well matched, albeit 
with predicted B aimer-Pickering absorption components that are 
too strong. Fig. [6] presents near infra-red spatially-resolved spec- 
troscopy of R136 WN stars (Schnurr et al. 2009) together with syn- 
thetic spectra, allowing for instrumental broadening (100 kms~^) 
plus additional rotational broadening of 200 kms~^ for R136a2, 
a3, c. As for NGC 3603 stars, the ratio of Hell 2.189/im to Br7 
provides an excellent diagnostic for the hydrogen content for the 
R136 stars. A summary of the resulting physical and chemical pa- 
rameters is presented in Table [S] with errors once again accounting 
for both photometric and spectroscopic uncertainties, representing 
the range of permitted values, although the primary uncertainty in- 
volves the distance to the LMC. 

With respect to earlier studies (Heap et al. 1994; de Koter 
et al. 1997; Crowther & Dessart 1998), the improved allowance 
for metal line blanketing also infers 20% higher stellar tempera- 
tures (T* ^53,000 K) for these stars, and in turn, larger bolomet- 
ric corrections (MboI — My ^ -4.6 mag). Such differences, with 
respect to non-blanketed analyses, are typical of Wolf-Rayet stars 
(Crowther 2007). Similar temperatures (T* = 50,000 or 56,000 K) 
to the present study were obtained for these WN stars by Riihling 
(2008) using a grid calculated from the Potsdam line-blanketed at- 
mospheric code (for a summary see Ruehling et al. 2008), albeit 
with 0.1-0.2 dex lower luminosities (logL/L© = 6.4 - 6.7) us- 
ing solely optically-derived (lower) extinctions and absolute mag- 
nitudes. 

To illustrate the sensitivity of mass upon luminosity, we use 
the example of R136a5 (02-3 IfAVN) for which we estimate T* -- 
50,000 K using the same diagnostics as the WN5h stars, in good 
agreement with recent analyses of 02 stars by Walborn et al. (2004) 
and Evans et al. (2010b). This reveals a luminosity of log L/Lq = 
6.35, corresponding to an initial mass in excess of 100 M©, versus 
T* r^ 42,500 K, log L/Lq = 5.9 and -65 M© according to de 
Koter et al. (1997). 

Regarding elemental abundances, we use the H/3 to He II 5412 
and Br7 to Hell 2.189/im ratios to derive (consistent) hydrogen 
contents, with the latter serving as the primary diagnostic for con- 



sistency with NGC 3603 stars (current observations do not include 
Hell 5412). We adopt scaled solar abundances for all metals other 
than CNO elements. Nitrogen abundances of Xat = 0.35% by mass, 
as predicted by evolutionary models, are consistent with near-IR 
N V 2.100/im recombination line observations, while we adopt car- 
bon and oxygen abundances of Xc = Xo = 0.004% by mass. 



4 EVOLUTIONARY MODELS 

We have calculated a grid of main- sequence models using the lat- 
est version of the Geneva stellar evolution code for 85, 120, 150, 
200, 300 and 500 M©. The main- sequence evolution of such high- 
mass stars does not suffer from stability issues. Although a detailed 
description is provided elsewhere (Hirschi et al. 2004), together 
with recent updates (Eggenberger et al. 2008), models include the 
physics of rotation and mass loss, which are both crucial to model 
the evolution of very massive stars. Although many details relating 
to the evolution of very massive stars remain uncertain, we solely 
consider the main- sequence evolution here using standard theoret- 
ical mass-loss prescription for O stars (Vink et al. 2001) for which 
Mokiem et al. (2007) provide supporting empirical evidence. In the 
models, we consider the onset of the Wolf-Rayet phase to take place 
when the surface hydrogen content Xh < 30% if Tes > 10,000 K, 
during which an empirical mass-loss calibration is followed (Nugis 
& Lamers 2000). The post-main sequence evolution is beyond the 
scope of this study and will be discussed elsewhere. 

All the main effects of rotation are included in the calcula- 
tions: centrifugal support, mass-loss enhancement and especially 
mixing in radiative zones (Maeder 2009), although predictions for 
both non-rotating and rotating models are considered here. For ro- 
tating models we choose an initial ratio of the velocity to critical 
(maximum) rotation of vinit/^^crit = 0.4, which corresponds to sur- 
face equatorial velocities of around 350 km s~^ for the 85 M© 
model and around 450 km s~^ for the 500 M© case. 

We have calculated both solar (Z=1.4% by mass) and LMC 
(Z=0.6% by mass) metallicities. The evolution of 120 M© models 
in the Hertzsprung-Russell diagram is presented in Fig[71 The ef- 
fects of rotation are significant. Due to additional mixing, helium 
is mixed out of the core and thus the opacity in the outer lay- 
ers decreases. This allows rotating stars to stay much hotter than 
non-rotating stars. Indeed, the effective temperature of the rotating 
models stay as high as 45,000 - 55,000 K, whereas the effective 
temperature of non-rotating models decreases to 20,000 - 25,000 
K. Therefore, rapidly rotating stars progress directly to the classi- 
cal Wolf-Rayet phase, while slow rotators are expected to become 
rj Car-like Luminous Blue Variables (see Meynet & Maeder 2005). 

Rotating models can reach higher luminosities, especially at 
very low metallicity (see Langer et al. 2007). This is explained by 
additional mixing above the convective core. Finally, by comparing 
models at solar metallicity and LMC metallicity, we can see that 
lower metallicity models reach higher luminosities. This is due to 
weaker mass-loss at lower metallicity. 

Langer et al. (2007) provide 150 M© tracks at 0.2 and 0.05 
Z© with high initial rotation velocities of 500 km s~^. From these 
it is apparent that rapidly rotating very metal-deficient models can 
achieve high stellar luminosities. According to Langer et al. (2007) 
a 150 M© model reaches log L/L© — 6.75 at 0.05 Z© and 6.5 
at 0.2 Z©. Could the R136 stars represent rapidly rotating stars 
of initial mass 150 M©? We have calculated a SMC metallicity 
(0.14 Z©) model for a 150 M© star initially rotating at i;init = 
450 kms~^, which achieves log L/L© ^ 6.5 after 1.5 Myr and 
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Figure 7. Comparison between main se- 
quence evolutionary predictions for both 
rotating (Z = 0.006: H120z06S400) and 
non-rotating (Z = 0.006: H120z06S000, Z 
= 0.014: H120zl4S000) 120 M© models. 
Horizontal lines in the top left-hand panel 
correspond to the transition to the Wolf- 
Rayet phase. 



6.67 after 2.5 Myr. Our models are thus compatible with Langer et 
al. (2007) described above. We nevertheless exclude the possibil- 
ity that the R136 WN5 stars possess initial masses below 150 M© 
since: 

(i) the metallicity of 30 Doradus is a factor of three higher than 
the SMC (e.g. Peimbert 2003, Lebouteiller et al. 2008); 

(ii) R136 has an age of less than 2 Myr (de Koter et al. 1998, 
Massey & Hunter 1998), since its massive stellar population is anal- 
ogous to the young star clusters in Car OBI (1-2 Myr, Walborn 
2010) and NGC 3603 (this study); 

(iii) Clumping-corrected mass-loss rates for the R136 WN stars 
agree well with LMC metallicity predictions (Vink et al. 2001), 
which are also supported by studies of O stars in the Milky Way, 
LMC and SMC (Mokiem et al. 2007). 

The surface abundances of a majority of high mass stars can 
be well reproduced by models of single stars including the effects 
of rotation. However, the VLT FLAMES survey has highlighted 
some discrepancies between models and observations (Hunter et 
al. 2008), questioning the efficiency of rotation induced mixing. We 
are currently investigating this matter by comparing the composi- 
tion of light elements like boron and nitrogen between models and 
observations and we find that models including rotation induced 
mixing reproduces the abundances of most stars well (Frischknecht 
et al. 2010). 

We can nevertheless consider how a less efficient rotation- 
induced mixing (or absence of mixing) would affect the conclu- 



sions of this paper. If rotation induced mixing was less efficient 
than our models predict, the masses derived here would remain at 
least as high, and would usually be higher. Indeed, less efficient 
mixing prevents the luminosity to increase as much since less he- 
lium is mixed up to the surface and the mean molecular weight, 
/i, remains lower (L ^ /nM^'^). In particular, less efficient mix- 
ing would prevent stars with initial masses around or below 150 
M0 from reaching such high luminosities as predicted in Langer 
et al. (2007, see discussion in previous paragraph) and would bring 
further support to our claim that these stars are more massive than 
150 M0. The challenge would then be to explain the high effective 
temperature (Tgff ^ 50 kK) observed for the R136 stars, which is 
best explained by rotation induced mixing. Note that less efficient 
mixing would also make impossible the quasi-chemical evolution 
of fast rotating (single or binary) stars (Yoon et al. 2006) which is 
currently one of the best scenarios for long/soft Gamma Ray Bursts 
progenitors. 



4.1 NGC 3603 

Fig.[8]compares various observational properties of NGC 3603 WN 
stars with solar metallicity evolutionary predictions. Non-rotating 



models imply current masses of \20Zir ^0 ^^^ 921^5 M© for 
Ala and Alb, respectively, at an age of ^^ 1.5 ± 0.1 Myr. These 
are in excellent agreement with dynamical mass determinations of 
116 ± 31 M0 and 89 =b 16 Mq for the primary and secondary Al 
components (Schnurr et al. 2008a). Initial and current stellar mass 
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Figure 8. Comparison between Solar metallicity (Z = 1.4%) models calculated for the main-sequence evolution of 85 - 200 Mq stars (initially rotating at 
Mnit/^crit =0-4 [dotted] and [solid]), and the physical properties derived from spectroscopic analysis of NGC 3603 WN6h stars. We obtain excellent 
agreement with dynamical masses for Ala and Alb for initially non-rotating models at ages of 1.5 ± 0.1 Myr. Current mass-loss rates match solar-metaUicity 
theoretical predictions (Vink et al. 2001) to within 0.2 dex. 
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Figure 9. Comparison between LMC-metallicity models calculated for the main-sequence evolution of 85 - 500 Mq stars, initially rotating at '^^init/'^crit = 
0.4 (dotted) or (solid) and the physical properties derived from our spectroscopic analysis. We obtain excellent agreement for initially rapidly rotating, 165 
- 320 M0 stars at ages of ~1.7 ± 0.2 Myr. Current mass-loss rates match LMC-metallicity theoretical predictions (Vink et al. 2001) to within 0.2 dex. 



estimates are shown in TableJH and include a (high) initial mass of 
166 zb 20 Mq for NGC 3603 B. Independent age estimates using 
pre-main sequence isochrones of low mass stars also favour low 
1±1 Myr ages (Sung & Bessell 2004), while Crowther et al. (2006) 
estimated 1.3±0.3 Myr for NGC 3603 from a comparison between 



massive O stars and theoretical isochrones (Lejeune & Schaerer 
2001). 
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Figure 10. Spectral comparison between VLT/SINFONI spectroscopy of 
N V 2.10/im in R136 WN5 stars (Schnurr et al. 2009) and synthetic spectra 
(red), allowing for instrumental broadening plus rotational broadening of 
200 kms-i for R136a2, a3 and c. 



4.2 R136 

Fig. |9] compares the derived properties of R136al, a2, a3 and c 
with LMC metallicity evolutionary predictions, under the assump- 
tion that these stars are single. Initial stellar masses in the range 
165 - 320 Mq are implied, at ages of 1.7 di 0.2 Myr, plus high 
initial rotational rates, in order that the observed surface hydrogen 
contents of 30 - 40% by mass are reproduced. Initial and current 
stellar mass estimates are included in Table \5\ Differences in age 
estimates reflect variations in initial rotation rates. Nevertheless, 
equatorial rotation rates of Ve ^ 200 (300) km s~^ are predicted 
after -1.75 Myr (2.75 Myr) for a 300 (150) M© star. 

In the absence of photospheric absorption features, the N V 
2.100/im feature provides the best diagnostic of rotational broaden- 
ing for the R136 stars. This recombination line is intrinsically nar- 
row at the SINFONI spectral resolution (FWHM — 15 A) because 
it is formed extremely deep in the stellar wind. Indeed, FWHM ^ 
15A is observed for R136al, as expected either for a non-rotating 
star, or one viewed pole-on (see Fig.[TO]). In contrast, a2 and a3 re- 
veal FWHM ^ 40A, corresponding to Ve sini ^ 200 kms~^. We 
are unable to quantify FWHM for R136c since the Nv feature is 
very weak, although it too is consistent with a large rotation rate. 
Therefore, a2, a3 and probably c show spectroscopic evidence for 
rapid rotation as shown in Fig.fTO] while al could either be a very 
slow rotator, or a rapid rotator viewed close to pole-on. 

For comparison, Wolff et al. (2008) found that R136 lacks 
slow rotators among lower mass 6-30 Mq stars. Wolff et al. de- 
rived Ve sin i = 189 ± 23 km s~^ for eleven 15-30 Mq stars within 
R136 versus ^;e sin i = 129 =b 13 km s~^ from equivalent mass field 
stars within the LMC. A much more extensive study of rotational 
velocities for O stars in 30 Doradus will be provided by the VLT- 
FLAMES Tarantula Survey (Evans et al. 2010a). 

Finally, although we defer the possibility that the WN5 stars 
in R136 are (equal mass) binaries until Sect. 5, it is necessary to re- 
mark upon the possibility of chance superpositions within the ob- 
served cluster. Line-of- sight effects should not play a significant 
role in our interpretation of bright systems such as R136al (Maiz 
Apellaniz 2008). Chance superposition of other stars has been cal- 
culated to contribute at most —10 to 20% of the (visible) light (J. 
Maiz Apellaniz, priv. comm.). For example, a 0.2 mag decrease in 
the absolute K-band magnitude of R136al arising from the contri- 
bution of lower mass stars along this sightline would lead to a 10% 



reduction in its initial (current) stellar mass, i.e. to 285 Mq (235 

Mo). 

4.3 Synthetic spectra 

We have calculated synthetic spectra for each of the LMC metal- 
licity models, both at the Zero Age Main Sequence (ZAMS) and 
ages corresponding to the surface hydrogen compositions reaching 
Xh = 30% for the rotating models. These are presented in Fig.fTTlin 
which mass-loss rates follow the theoretical mass-loss recipes from 
Vink et al. (2001), both for the case of radially-dependent clumped 
winds with a volume filling factor of /oo=0.1 at ^;oo (solid) and 
smooth (/oo=l, dotted) winds. For clumped winds, these result in 
ZAMS synthetic O supergiant spectra which are equivalent to 03 If 
(for 85 Mo) and 02-3IfAVN5-6 (for 200 Mo) subtypes. Weaker 
He II A4686 emission would naturally be predicted if we were to 
instead adopt smooth winds, including 03 V for the 85 Mo case 
and 02 Ill-If for the 200 Mo case at the zero age main sequence. 

Although specific details depend upon the degree of wind 
clumping in O star winds and validity of the smooth Vink et al. 
predictions, stars whose masses exceed a given threshold are ex- 
pected to display a supergiant signature from the outset. Overall, 
these results suggest a paradigm shift in our understanding of early 
O dwarfs, which have hitherto been considered to represent ZAMS 
stars for the highest mass stars (Walborn et al. 2002). 

For the NGC 3603 and R136 WN stars studied here, wind 
clumping is required to reproduce the electron scattering wings of 
He II A4686 (Hillier 1991). Overall, we find log M - log Mvink = 
+0A3 ± 0.09, in close agreement with Martins et al. (2008) who 
found that spectroscopic mass-loss rates for the Arches WN7-9 
stars exceeded predictions by -1-0.2 dex. 

Very powerful winds naturally result from the dependence of 
the mass-loss rate upon the ratio of radiation pressure to gravity, 
Ve oc L/M. Since L oc M^'^ for ZAMS stars in the range 30 
- 300 Mo based upon evolutionary calculations discussed above, 
Fe oc M°^. A 30 Mo main-sequence star possesses an Edding- 
ton parameter of Fe -- 0.12 (Conti et al. 2008), therefore a 300 
Mo star will possess Fe ^ 0.5. The latter achieve Fe ^ 0.7 after 
1.5 Myr, typical of the cluster age under investigation here. Fe in- 
creases from 0.4 to 0.55 for stars of initial mass 150 Mo after 1.5 
Myr. 

Therefore, the very highest mass stars in very young systems 
may never exhibit a normal O-type absorption line spectrum. In- 
deed, 02-3 V stars may rather be limited to somewhat lower ZAMS 
stars, with a higher mass threshold at lower metallicity. Recall the 
relatively low dynamical mass of 57 Mo for the 03 V primary in 
Rl 36-38 (Massey et al. 2002), which currently represents the most 
massive 02-3 V star dynamically weighed, versus 83 Mo (WR20a: 
Bonanos et al. 2004; Rauw et al. 2005), >87 Mo (WR21a: Niemela 
et al. 2008) and 1 16 Mo (NGC 3603 Ala: Schnurr et al. 2008a) for 
some of the most massive H-rich WN stars. O If* or even WN-type 
emission line spectra may be expected for the highest mass stars 
within the very youngest Giant H II regions, as is the case for R136 
and NGC 3603 here, plus Car OBI (Smith 2006) and W43 (Blum 
et al. 1999). 



5 BINARITY 

NGC 3603 Al and C are confirmed binaries (Schnurr et al. 2008a), 
and R136c is a probable binary (Schnurr et al. 2009) while compo- 
nent B and R136al, a2, a3 are presumed to be single. Of course. 
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Figure 11. Synthetic visual spectra calculated for zero-age main-sequence stars (black) of initial mass 85 and 200 Mq, initially rotating at '^^init/'^crit = 
or 0.4, assuming that winds are clumped (/oo =0.1, solid lines) and follow Vink et al. (2001). Synthetic spectra are also shown at the time at which Xh = 
30% by mass is predicted at the surface of the initially rotating model, i.e. 1.75 - 2.75 Myr (red). Synthetic spectra are broadened by 50 kms"-*^ ('^init /'^^crit 
= 0), Ve sini = 250 kms~^ ('i^init/'^^crit = 0.4, 0.01 Myr), Ve sini = 150 kms~^ ('^init/'^crit = 0.4, 1.75-2.75 Myr). Heii A4686 emission is predicted for 
all phases for these clumpy winds - most prominently during the later phases of initially rapidly rotating very massive stars. If we were to adopt non-clumped 
winds (/oo = 1, dotted lines), Heii A4686 absorption would be predicted for the ZAMS 85 Mq models. 



we cannot unambiguously confirm that these are genuinely single 
but if their mass ratios differ greatly from unity the derived physi- 
cal properties reflect those of the primary. Therefore, here we shall 
focus upon the possibility that the R136 stars are massive binaries, 
whose ratio is close to unity (Moffat 2008). 



If we adopt similar mean molecular weights, /i, for individual 
components, current (evolutionary) masses of 150 -i- 150 Mq, 200 
-1-100 Mq and 220 -i- 55 Mq would be required to match the cur- 
rent properties of R136al for mass ratios of 1, 0.5 and 0.25 since 
L oc /iM^^. Of these possibilities, only near-equal mass binaries 
would contradict the high stellar masses inferred here. If the sepa- 
ration between putative binary components were small, radial ve- 
locity variations would be expected. There is no unambiguous ev- 
idence for binarity among the stars under investigation, although 
R136c does exhibit marginal radial velocity variability (Schnurr et 
al. 2009). 



In contrast, the WN6h system Al in the young Milky Way 
cluster NGC 3603 is a short-period massive binary (Schnurr et al. 
2008a). However, longer period systems with periods of months to 
years could have easily avoided spectroscopic detection. We shall 
employ anticipated properties of colliding wind systems and dy- 
namical interactions to assess whether the R136a stars are realisti- 
cally long-period binaries. This approach is especially sensitive to 
equal wind momenta (equal mass) systems since the X-ray emis- 
sion is maximal for the case of winds whose momenta are equal. Al 
in NGC 3603 is relatively faint in X-rays because the components 
are in such close proximity that their winds collide at significantly 
below maximum velocity. 



Table 6. Predicted X-ray luminosities from R136al for scenarios in which 
it is either a colliding wind binary using analytical predictions from Stevens 
et al. (1992) and Pittard & Stevens (2002), or single/multiple for various 
empirical Lx/i^Boi values. Analytical calculations adopt component sepa- 
rations of d = 3, 30, 300 AU, for mass ratios: (i) g = 1 (150 + 150 Mq) and 
equal wind momenta, 77 = 1 (Mi = 2.8 x 10~^ Mq yr~i, -Uoo,! = 2600 
kms-i); (ii) q = 0.25 (220 + 55 Mq) and ry = 0.05 (Mi = 4.2 x IQ-^ 
MQyr-i,i;oo,i=2600kms-i,M2 = 3x IQ-^ Mq yr-^, ^00,2 =2000 
kms~^). The intrinsic X-ray luminosity of R136a is taken from Guerrero 
& Chu (2008). 
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5.1 X-rays 

If both components in a binary system possessed similar masses, 
once outflow velocities of their dense stellar winds achieve asymp- 
totic values, collisions would produce stronger X-ray emission 
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than would be expected from a single star. Analytical estimates 
of X-ray emission from colliding wind systems are approximate 
(Stevens et al. 1992). Nevertheless, this approach does provide con- 
straints upon orbital separations to which other techniques are cur- 
rently insensitive and is especially sensitive to systems whose wind 
strengths are equal. According to Pittard & Stevens (2002) up to 
17% of the wind power of the primary (and 17% of the secondary) 
can be radiated in X-rays for equal wind momenta systems, ver- 
sus 0.4% of the primary wind power (56% of the secondary) for 
(unequal mass) systems whose wind momentum ratio is 77 = 0.01. 

To illustrate the diagnostic potential for X-ray observations, 
let us consider the expected X-ray emission from R136c under 
the assumption that it is single. The intrinsic X-ray luminosities 
of single stars can be approximated by Lx/LboI ^^ 10 ~^ (Chle- 
bowski et al. 1989). From our spectroscopic analysis we would 
expect Lx ~ 2 X 10^^ ergs~^, yet Chandra imaging reveals an 
intrinsic X-ray luminosity which is a factor of 30-50 times higher 
(Portegies Zwart et al. 2002; Townsley et al. 2006; Guerrero & Chu 
2008), arguing for a colliding wind system in this instance. Ana- 
lytical models favour equal mass components separated by '-^lOO 
Astronomical Units (AU). 

In contrast, the expected X-ray emission from the sum of 
R136al, a2 a3 (and a5) - unresolved at Chandra resolution - is 
Lx ~ 7 X 10^^ erg s~^ under the assumption that they are single. 
In this case, the intrinsic X-ray emission from R136a is observed to 
be only a factor of ^3 higher (Portegies Zwart et al. 2002; Guerrero 
& Chu 2008) arguing against short period colliding wind systems 
from any R136a components. Multiple wind interactions (outside 
of the binaries but within the cluster) will also produce shocks and 
X-ray emission (see e.g. Reyes-Iturbide et al. 2009). 

Empirically, colliding winds within 0-type binaries typically 
exhibit Lx/LboI "^ 10~^ (Rauw et al. 2002), although binaries 
comprising stars with more powerful winds (i.e. Wolf-Rayet stars) 
often possess stronger X-ray emission. Indeed, NGC 3603 C has an 
X-ray luminosity of Lx ^ 4 x 10^^ ergs"^ (Moffat et al. 2002), 
corresponding to Lx/LboI ^ 5 x 10~^ based upon our spectroscop- 
ically derived luminosity (similar results are obtained for R 136c). If 
we were to assume this ratio for the brightest components of R136a 
(al, a2, a3 and a5) under the assumption that each were close bina- 
ries - i.e. disregarding predictions from colliding wind theory - we 
would expect an X-ray luminosity that is a factor of 15 times higher 
than the observed value. Therefore, the WN stars in R136a appear 
to be single, possess relatively low-mass companions or have wide 
separations. 

If we make the reasonable assumption that equal mass stars 
would possess similar wind properties, we can consider the effect 
of wind collisions upon the production of X-rays. Let us consider 
a 150 -I- 150 M© binary system with a period of 100 days in a 
circular orbit with separation 3 AU, whose individual components 
each possess mass-loss rates of M = 2.8 x 10~^ M© yr~^ and 
wind velocities of ^;oo = 2600 km s~^. A pair of stars with such 
properties could match the appearance of R136al, with individual 
properties fairly representative of R136a3. 

Since the ratio of their wind momenta, 77 is unity, the fraction, 
S, of the wind kinetic power processed in the shock from each star 
is maximal (E ^1/6, Pittard & Stevens 2002). One must calculate 
the conversion efficiency, x^ of kinetic wind power into radiation 
for each star. From Stevens et al. (1992) 



X' 



vtdi2 

M_7 ' 



where vg is the wind velocity in units of 10^ cms ^, di2 is the 



distance to the interaction region in units of 10^^ cm (i.e. half the 
separation for equal wind momenta) and M-7 is the mass-loss rate 
in units of \0~^ Mq yr~^ . We set a lower limit of x = 1 for cases in 
which the system radiates all of the collision energy. This allows an 
estimate of the intrinsic X-ray luminosity (J. Pittard, priv. comm.). 

For a 150 -I- 150 M© system with separation of 3 AU, we derive x "^ 
3.7 and predict Lx = 2.7 x 10^^ erg s~^ for each component. The 
total exceeds the measured X-ray luminosity of R136a by a factor 
of -200 (Townsley et al. 2006; Guerrero & Chu 2008), as shown 
in Table [6l Recall that if R136al were to be a massive binary with 
a mass ratio of 0.25, the secondary would possess a mass of —55 
M0. This would contribute less than 10% of the observed spectrum, 
with a secondary to primary wind momentum ratio of 77 —0.05. 
Predictions for this scenario are also included in Table [6] revealing 
X-ray luminosities a factor of —6 lower than the equal wind (equal 
mass) case. 

Stevens et al. (1992) also provide an expression for the char- 
acteristic intrinsic column density in colliding wind binaries, A^h 
(their Eqn. 11), from which Nu — 2.4 xlO^^ cm~^ would be in- 
ferred for R136al is it were an equal mass binary system, separated 
by 3 AU. This is an order of magnitude higher than the estimate for 
R136a (Townsley et al. 2006; Guerrero & Chu 2008), arising from 
a combination of foreground and internal components. 

Of course, a significant fraction of the shock energy could pro- 
duce relativistic particles rather than X-ray emission (e.g. Pittard & 
Dougherty 2006). Formally, an equal mass binary system whose 
components are separated by 600 AU is predicted to produce an 
X-ray luminosity of Lx = 2.7 x lO'^'^ erg s~^, which is comparable 
to the total intrinsic X-ray luminosity of R136a (Guerrero & Chu 
2008). If 90% of the shock energy were to contribute to relativis- 
tic particles, the same X-ray luminosity could arise from a system 
whose components were separated by 60 AU. Still, the relatively 
low X-ray luminosity of R 136a suggests that if any of the R136a 
WN5 stars are composed of equal mass binaries, their separations 
would need to be in excess of 100-200 AU. How wide could very 
massive binaries be within such a dense cluster? To address this 
question we now consider dynamical interactions. 



5.2 Dynamical interactions 

The binding energy of a 'hard' 100 Mq -i- 100 Mq binary with a 
separation of 100 AU is comparable to the binding energy of a mas- 
sive cluster ( 10^^ J). Such a system will have frequent encounters 
with other massive stars or binaries. The collision rate for a sys- 
tem of separation, a, in a cluster of number density, n, and velocity 
dispersion, a, is 

Tcoii -30ncra^(l - 0) 

where is the Safranov number, indicating the importance of grav- 
itational focusing (Binney & Tremaine 1987). Obviously, an en- 
counter with a low mass star will not affect a very massive binary, 
so we assume only that encounters with stars in excess of 50 Mq 
will significantly affect a very massive binary. 

For a 100 -I- 100 M© binary with 100 AU separation in a cluster 
with a velocity dispersion of a = 5 km s~^, initially containing 
a number density of 500 pc"^ (n = 1.5 x 10"'^'^ m"^) then = 
35. Each very massive binary should have an encounter every 1.8 
Myr. For binaries wider than 100 AU the encounter rate will be 
even higher, since the encounter rate scales with the square of the 
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Table 7. Comparison between the properties of the NGC 3603, Arches and 
R136 star clusters and their most massive stars, Mmrx- 



Name 



Mel r 
Mq Myr 



Ref 



lyrinit 

Mq 



/\/fcurrent 

Mq 



Ref 



NGC 3603 lO'^ -1.5 a, b 166±20 132±13 b 

Arches 2x10^ 2.5±0.5 c, d 120-150 >95-120 b, d 

~1.8 b >185l^^ >130l25 t> 

R136 <5.5xl0'^ -1.7 e, b 320+J°° 265+^° b 



(a) Harayama et al. (2008); (b) This work; (c) Figer (2008); (d) Martins et 
al. (2008); (e) Hunter et al. (1995) 

separation, i.e. encounters every ^0.2 Myr for separations of 300 
AU or 0.05 Myr for separations of 600 AU. 

Such close encounters will create an unstable multiple system 
that will rapidly decay by ejection of the lowest mass star (Anosova 
1986). This would not necessarily destroy the binary. Instead, it will 
harden the system in order to gain the energy required to eject the 
other star. Thus very massive systems separated by more than '-^lOO 
AU will reduce their separations through dynamical interactions 
with other high mass stars within the dense core of R136a. 

Should one of the high mass binaries in R136a elude harden- 
ing in the way outlined above, would it still escape detection if it 
remained at a large separation? One binary with components of 150 
+ 150 Mq separated by -300 AU for which -30% of the shock 
energy contributes to the X-ray luminosity, together with intrinsic 
X-ray luminosities from single R136a members, obeying Lx/LboI 
— 10~^ (Chlebowski et al. 1989) would indeed mimic the observed 
R136a X-ray luminosity (Guerrero & Chu 2008). 

In summary, if we adopt similar ratios of X-ray to bolomet- 
ric luminosities as for R136c and NGC 3603C, we would expect 
an X-ray luminosity from R136a that is a factor of 15 times higher 
than the observed value if al, a2 and a3 were each colliding wind 
systems. Alternatively, we have followed the colliding wind theory 
of Stevens et al. (1992) and Pittard & Stevens (2002), and conser- 
vatively assume 30% of the shock energy contributes to the X-ray 
luminosity. If all the very massive stars within the Chandra field-of- 
view (R136al, a2, a3 and a5) were members of equal mass binaries 
with separations of 300 AU, we would expect a X-ray luminosity 
that is over a factor of 2 higher than observed. In view of dynami- 
cal effects, at most one of the WN components of R136a might be 
a long period, large separation (>300 AU) equal mass binary. We 
cannot rule out short-period, highly unequal-mass binary systems 
of course, but such cases would have little bearing upon our derived 
stellar masses. 



6 CLUSTER SIMULATIONS 

We now compare the highest mass stars in R136 and NGC 3603 
with cluster predictions, randomly sampled from the stellar IMF, 
for a range of stellar mass limits. For R136, the stellar mass within 
a radius —4.7 pc from R136al has been inferred for >2.8 M© at 
2.0 X 10^ Mq (Hunter et al. 1995). If the standard IMF (Kroupa 
2002) is adopted for < 2.8 M© a total cluster mass of <5.5 x 10^ 
Mq is inferred. Andersen et al. (2009) estimated a factor of two 
higher mass by extrapolating their measured mass down to 2. 1 Mq 
within 7 pc with a Salpeter slope to 0.5 Mq. These are upper limits 
to the mass of R136 itself since the mass function in the very cen- 
tre of R136 cannot be measured (Maiz Apellaniz 2008), such that 
they include the associated halo of older (> 3 Myr) early- type stars 



viewed in projection - indeed Andersen et al. (2009) obtained an 
age of 3 Myr for this larger region. Examples of such stars include 
R136b (WN9h), R134 (WN6(h)) and Melnick 33Sb (WC5). 

For NGC 3603, we adopt a cluster mass of 1.0 x 10^ M©, the 
lower limit obtained by Harayama et al. (2008) from high resolu- 
tion VLT NAOS/CONICA imaging. These are presented in Table|7l 
together with the Arches cluster, for which Figer (2008) estimated 
a cluster mass of 2 x 10^ Mq based on the mass function of Kim 
et al. (2006). 



6.1 Cluster populations 

We simulate a population of clusters and stars by randomly sam- 
pling first from a power-law cluster mass function (CMF), and then 
populating each cluster with stars drawn randomly from a stellar 
IMF (Parker & Goodwin 2007). Cluster masses are selected from 
a CMF of the form N{M) oc M~^, with standard slope /3 = 2 
(Lada & Lada 2003) between cluster mass limits 50 Mq and 2 x 
10^ M0. These limits enable the full range of cluster masses to be 
sampled, including those with masses similar to that of R136. The 
total mass of clusters is set to 10^ Mq to fully sample the range of 
cluster masses. Each cluster is populated with stars drawn from a 
three-part IMF (Kroupa 2002) of the form 



N{M) 



oc 



M+°-3 mo < M/Mq < mi , 
M"^-^ mi < M/Mq <m2, 
■^•^ m2 < M/Mq < ms , 



M 



where mo = 0.02 M©, mi = 0.1 M© and m2 = 0.5 M©. 
We use three different values for ms in the simulations; 150 Mq, 
300 Mq and 1000 Mq. Stellar mass is added to the cluster until the 
total mass is within 2% of the cluster mass. If the final star added 
to the cluster exceeds this tolerance, then the cluster is entirely re- 
populated (Goodwin & Pagel 2005). 

Our random sampling of the IMF allows low-mass clusters 
to be composed of one massive star and little other stellar material. 
This contravenes the proposed fundamental cluster mass-maximum 
stellar mass (CMMSM) relation (Weidner & Kroupa 2006, Weid- 
ner et al. 2010). However, the average maximum stellar mass for a 
given cluster mass closely follows the CMMSM relation (Parker & 
Goodwin 2007, Maschberger & Clarke 2008). The results for our 
three Monte Carlo runs are shown in Fig[T2] They reveal that the av- 
erage relation between the maximum stellar mass and cluster mass 
(in the 10^ and 10^ Mq interval) is recovered (Weidner et al. 2010), 
without the constraint that the cluster mass governs the maximum 
possible stellar mass (Parker & Goodwin 2007). We also indicate 
25% and 75% quartiles for the instances of a 10^ M© (NGC 3603), 
2 X 10^ Mo (Arches) and 5 x 10^ M© (R136) cluster. 



6.2 R136 

If we assume that R136al, a2 and c are either single or the primary 
dominates the optical/IR light, we obtain an average of 260 Mq for 
their initial masses. In reality this value will be an upper limit due to 
binarity and/or line-of- sight effects. Nevertheless, from Fig.[T2l we 
would expect the average of the three most massive stars in a cluster 
of mass 5 x 10^ M© to be 150 Mq, 230, 500 M© for an adopted 
upper mass limit of ms = 150, 300 and 1000 M©, respectively. 
Therefore, an upper limit close to 300 M© is reasonably consistent 
with the stellar masses derived (see also Oey & Clarke 2005). 

In the lower panel of Figure [13] we therefore present typical 
mass functions for a cluster of mass 5x10^ M© for an adopted 
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Figure 12. Average mass of the highest mass star (solid black line) and 25% and 75% quartiles for 10^, 2 x 10^ and 5 x 10^ Mq (vertical error bars), plus the 
average of the three highest mass stars (dotted red line) versus cluster mass. Clusters were randomly populated by sampling a three-part IMF (Kroupa 2002) 
using a Monte Carlo simulation, under the assumption that the upper mass limit is (left) 150 Mq; (middle) 300 Mq; (right) 1000 Mq. In order to minimise 
sampling problems we take an 'average' cluster as the mean of 1000 realisations of a cluster of a particular mass. 



upper limit of ms = 300 Mq. From a total of 1.05-1.10 x 10^ 
stars we would expect ^14 initially more massive than 100 M0 
to have formed within 5 parsec of R136al, since this relates to the 
radius used by Hunter et al. (1995) to derive the cluster mass. In- 
deed, of the 12 brightest near-infrared sources within 5 parsec of 
R136al, we infer initial evolutionary masses in excess of 100 Mq 
for 10 cases (all entries in Table [5] except for R136b and R139). Of 
course, this comparison should be tempered by (a) contributions of 
putative secondaries to the light of the primary in case of binarity; 
(b) dynamical ejection during the formation process (e.g. Brandl et 
al. 2007). Indeed, Evans et al. (2010b) propose that an 02III-If* 
star from the Tarantula survey is a potential high velocity runaway 
from R136, in spite of its high stellar mass (^ 90Mq). 

Overall, R136 favours a factor of ^2 higher stellar mass limit 
than is currently accepted. We now turn to NGC 3603 to assess 
whether it supports or contradicts this result. 

6.3 NGC 3603 

On average, the highest initial stellar mass expected in a 10^ Mq 
cluster would be ^ 120 M© with an upper mass limit of ms = 150 
Mo, - 200 Mq for ms = 300 M© and - 400 M© for ms = 10^ 
M0. Since we infer an initial mass of 1 66^20 ^0 for component 
B, a stellar limit intermediate between 150 and 300 Mq might be 
expected, if it is either single or dominated by a single componeno 
However, no stars initially more massive than 150 Mq would be 
expected in 25% of 10^ Mq clusters for the case of ms = 300 Mq 
(recall middle panel of Fig.[T2]). 

In the upper panel of Figure [13] we present typical mass func- 
tions for a cluster of mass 1x10^ Mq for an adopted upper mass 
limit of ms = 300 Mq. We would expect 3 stars initially more mas- 
sive than 100 M©, versus 4 observed in NGC 3603 (TableEJ. If we 
were to extend the upper limit of the IMF to ms = 1000 Mq, then 
the average of the three most massive stars of a 10^ Mq cluster 
would be 230 M© which is not supported by NGC 3603 (150 Mq). 

Therefore, the very massive stars inferred here for both R136 
and NGC 3603 are fully consistent with a mass limit close to ms 
= 300 Mq, both in terms of the number of stars initially exceeding 
100 Mq and the most massive star itself. Before we conclude this 
section, let us now consider the Arches cluster. 



■^ We do not attempt to claim that NGC 3603 B is single, even though it is 
photometrically and spectroscopically stable and is not X-ray bright (Moffat 
et al. 2004; Schnurr et al. 2008a). 



6.4 Arches 

The upper mass limit of ^150 M© from Figer (2005) was obtained 
for the Arches cluster from photometry. This was confirmed from 
spectroscopic analysis by Martins et al. (2008) who estimated ini- 
tial masses of 120 - 150 M© range for the most luminous (late WN) 
stars. How can these observations be reconciled with our results for 
NGC 3603 and R136? 

Let us consider each of the following scenarios: 

(i) the Arches cluster is sufficiently old (2.5 ± 0.5 Myr, Martins 
et al. 2008) that the highest mass stars have already undergone core- 
collapse, 

(ii) the Arches is unusually deficient in very massive stars for 
such a high mass cluster, or 

(iii) Previous studies have underestimated the masses of the 
highest mass stars in the Arches cluster. 

We have compared the derived properties for the two most 
luminous WN stars in the Arches cluster (F6 and F9) from Mar- 
tins et al. (2008) to the solar metallicity grids from Sect. 4. These 
are presented in Table [U revealing initial masses of > 115^25 
Mq from comparison with non-rotating models. Uncertainties 
from solely K-band spectroscopy are significantly higher than our 
(UV)/optical/K-band analysis of WN stars in NGC 3603 and R136. 
Masses represent lower limits since Martins et al. (2008) indicated 
that the metallicity of the Arches cluster is moderately super- solar. 

Core hydrogen exhaustion would be predicted to occur after 
^-^2.5 Myr for the non-rotating 150 M© solar metallicity models, 
with core-collapse SN anticipated prior to an age of 3 Myr. How- 
ever, comparison between the properties of F6 and F9 - the most 
luminous WN stars in the Arches cluster - with our solar metallic- 
ity evolutionary models suggest ages of ^2 Myr. On this basis, it 
would appear that scenario (i) is highly unlikely. 

Let us therefore turn to scenario (ii), on the basis that the initial 
mass of the most massive star within the Arches cluster was 150 
Mq. On average, the highest mass star in a 2 xlO^ M© cluster 
would be expected to possess an initial mass in excess of 200 M© 
for an upper mass limit of ms = 300 Mq. From our simulations, 
the highest mass star within such a cluster spans a fairly wide range 
(recall Fig. [12]). However, an upper limit of 120 (150) M© would 
be anticipated in only 1% (5%) of cases if we were to adopt ms = 
300 Mq . Could the Arches cluster be such a statistical oddity? 

Our simulations indicate that a 2 x 10^ M© cluster would be 
expected to host six stars initially more massive than 100 M© for 
ms = 300 Mq. Our solar metallicity evolutionary models com- 
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Table 8. Comparison between stellar properties (and inferred masses) of the most luminous Arches stars from Martins et al. (2008) obtained using different 
photometry, foreground interstellar extinctions and distances to the Galactic Centre. 



Name 


SpType 


niKs 


niH - niKs 


Ref 


Ak, 


Ref 


(m-M)o 


Ref 


Mk, 


BCk. 


Ref 


logL 


Minit 


Mcurrent 


Ref 






mag 


mag 




mag 




mag 




mag 


mag 




Lq 


Mo 


Mo 




F6 


WN8-9h 


10.37* 


1.68* 


a 


2.8 ±0.1 


b 


14.4 ±0.1 


c 


-6.83 ±0.14 


-4.0 ± 0.35 


d 


6.25 ±0.15 


>n5tll 


> 90±1^ 


h 






10.07 


1.96 


e 


3.1 ±0.2 


f 


14.5 ±0.1 


g 


-7.53 ± 0.22 


-4.0 ± 0.35 


d 


6.53 ±0.16 


>185±p3 


>130+^^ 


h 


F9 


WN8-9h 


10.77* 


1.67* 


a 


2.8 ±0.1 


b 


14.4 ±0.1 


c 


-6.43 ±0.14 


-4.4 ± 0.35 


d 


6.25 ±0.15 


>n5tll 


> 92t\l 


h 






10.62 


1.78 


e 


3.1 ±0.2 


f 


14.5 ±0.1 


g 


-6.98 ± 0.22 


-4.4 ± 0.35 


d 


6.47 ±0.16 


>165±fo 


>n()tll 


h 



(a) Figer et al. (2002); (b) Stolte et al. (2002); (c) Eisenhauer et al. (2005); (d) Martins et al. (2008); (e) Espinoza et al. (2009); (f) Kim et al. (2006); (g) Reid 

(1993); (h) This work 
J: NICMOS F205W and F160W filters; *: Includes -0.25 mag offset to bolometric correction relative to F. Martins (priv. comm.), as reported in Clark et al. 

(2009). 



bined with spectroscopic results from Martins et al. (2008) sug- 
gest that the 5 most luminous WN stars (Fl, F4, F6, F7, F9), with 
log L/Lq > 6.15, are consistent with initial masses of > 100 Mq, 
providing they are single or one component dominates their near-IR 
appearance. 

Finally, let us turn to scenario (iii), namely that the stellar 
masses of the Arches stars have been underestimated to date. Mar- 
tins et al. (2008) based their analysis upon HST/NICMOS photom- 
etry from Figer et al. (2002) together with a (low) Galactic Centre 
distance of 7.6 kpc from Eisenhauer et al. (2005) plus a (low) fore- 
ground extinction of Aks =2.8 mag from Stolte et al. (2002). Let 
us also consider the resulting stellar parameters and mass estimates 
on the basis of the standard Galactic Centre distance of 8 kpc (Reid 
1993), more recent (higher) foreground extinction of Ak^ = 3.1 ± 
0.19 mag from Kim et al. (2006)^lus VLT/NACO photometric 
results from Espinoza et al. (2009o These yield absolute K-band 
magnitudes that are 0.55-0.7 mag brighter than Figer et al. (2002), 
corresponding to 0.22-0.28 dex higher bolometric luminosities. 
With respect to the hitherto 150 Mo stellar mass limit identified by 
Figer (2005), these absolute magnitude revisions would conspire to 
increasing the limit to >200 Mo- 

In Table [8] we provide the inferred properties of F6 and F9 
on the basis of these photometric properties plus the stellar tem- 
peratures and K-band bolometric corrections derived by Martins et 
al. (2008). Our solar metallicity non-rotating models indicate ini- 
tial masses of > 185^1^ Mo and ^6^165^^° Mo for F6 and F9, 
respectively. These are lower limits to the actual initial masses, in 
view of the super- solar metallicity of the Arches cluster (Martins 
et al. 2008). In total, ^5 stars are consistent with initial masses in 
excess of :^150 Mo (those listed above), plus a further 5 stars for 
which initial masses exceed -100 Mo (F3, F8, F12, F14, F15). 

Overall, ten stars initially more massive than '-^lOO would sug- 
gest either that the mass of the Arches cluster approaches 3x10^ 
Mo or several cases are near-equal mass binaries (Lang et al. 2005, 
Wang et al. 2006). As such, we would no longer require that the 
Arches cluster is a statistical oddity, since the highest mass star 
would not be expected to exceed 200 Mo in 25% of 2 x 10^ Mo 
clusters. 

In conclusion, we have attempted to reconcile the properties of 
the highest mass stars in Arches cluster with an upper mass limit of 
7713 = 300 Mo- Based upon the Martins et al. (2008) study, com- 



2 Indeed, if we adopt (H-Ks)o = -0.1 1 mag (Crowther et al. 2006) for F6 
and F9, Axg = 3.2 ± 0.2 mag is implied from the Espinoza et al. (2009) 



photometry. 



parison with evolutionary models suggests that the initial masses of 
the most massive stars are > 115^25 ^0 with ages of —2 Myr. We 
find that an upper mass of 120 Mo would be expected in only 1% 
of 2 X 10^ Mo clusters for which ms = 300 Mo- However, use of 
contemporary near-IR photometry and foreground extinctions to- 
wards the Arches cluster, together with the standard 8 kpc Galactic 
Centre distance reveal an initial mass of > I85I45 ^0 ^^^ ^^e most 
massive star based on solar metallicity evolutionary models, with 
4-5 stars consistent with initial masses >150 Mo- An upper mass 
of 200 Mo would be expected in 25% of cases. Robust inferences 
probably await direct dynamical mass determinations of its bright- 
est members, should they be multiple. 



7 GLOBAL PROPERTIES OF R136 

We will re-evaluate the ionizing and mechanical wind power result- 
ing from all hot, luminous stars in R136 and NGC 3603 elsewhere 
(E. Doran et al. in preparation). Here, we shall consider the role 
played by the very massive WN stars in R136. We have updated the 
properties of early-type stars brighter than My = -4.5 mag within 
a radius of ^^5 parsec from R136al (Crowther & Dessart 1998) 
to take account of contemporary Tgff-spectral type calibrations for 
Galactic stars (Conti et al. 2008), and theoretical mass-loss rates 
(Vink et al. 2001) for OB stars. Until a census of the early- type 
stars within R 136 is complete, we adopt 03 subtypes for those stars 
lacking spectroscopy (Crowther & Dessart 1998). Aside from OB 
stars and the WN5h stars discussed here, we have included the con- 
tribution of other emission-line stars, namely 02-3 IfAVN stars for 
which we adopt identical temperatures to those of the WN5 stars, 
one other WN5h star (Melnick 34), a WN6(h) star (R134) for which 
we estimate T* ^ 42,000 K and a WN9h star (R136b) for which 
we estimate T* ^ 35,000 K from its HST/FOS and VLT/SINFONI 
spectroscopy. Finally, a WC5 star (Melnick 33Sb) lies at a projected 
distance of 2.9 pc for which we adopt similar wind properties and 
ionizing fluxes to single WC4 stars (Crowther et al. 2002). 

Fig.[T4]shows the integrated Lyman continuum ionizing fluxes 
and wind power from all early-type stars in R136, together with 
the explicit contribution of R136al, a2, a3 and c, amounting to, 
respectively, 46% and 35% of the cluster total. If we were to adopt 
a 2,000 K systematically higher temperature calibration for 0-type 
stars in the LMC, as suggested by some recent results (Mokiem et 
al. 2007, Massey et al. 2009), the contribution of WN5h stars to 
the integrated Lyman continuum flux and wind power is reduced to 
43% and 34%, respectively. 
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Figure 13. Theoretical mass functions for clusters with stellar masses of 
(upper) 10^ M0 and (lower) 5 xlO^ M©, for the scenario with an upper 
mass limit of ms = 300 Mq. 3 and 14 stars with initial masses in excess of 
100 M0 are anticipated, respectively. Our Monte Carlo approach (Parker & 
Goodwin 2007) breaks down for very massive stars in the upper panel due 
to small number statistics. 
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Figure 14. (Upper panel) Incremental Lyman continuum ionizing radi- 
ation from early-type stars within 5 parsec from R136al (Crowther & 
Dessart 1998) either following the Galactic O subtype-temperature calibra- 
tion (solid black), or systematically increasing the temperature calibration 
by 2,000 K for LMC stars (dotted line). The four very luminous WN5h stars 
discussed here (dashed red line) contribute 43 - 46% of the total; (Lower 
panel) As above, except for the mechanical wind power from early-type 
stars within 5 parsec of R136al, based upon theoretical wind prescriptions 
(Vink et al. 2001) to which the four WN5h stars contribute 34% of the total. 



We have compared our empirical results with population syn- 
thesis predictions (Leitherer et al. 1999) for a cluster of mass 5.5 
X 10^ M0 calculated using a standard IMF (Kroupa 2002) and 
evolutionary models up to a maximum limit of 120 Mq. These 
stars alone approach the mechanical power and ionizing flux pre- 
dicted for a R136-like cluster at an age of ^1.7 Myr (Leitherer et 
al. 1999). Indeed, R136al alone provides the Lyman continuum 
output of seventy 07 dwarf stars (Conti et al. 2008), supplying 
7% of the radio-derived '-^lO^^ photon s~^ ionizing flux from the 
entire 30 Doradus region (Mills et al. 1978, Israel & Koornneef 
1979). Improved agreement with synthesis models would be ex- 
pected if evolutionary models allowing for rotational mixing were 
used (Vazquez et al. 2007), together with contemporary mass-loss 
prescriptions for main- sequence stars. 



8 DISCUSSION AND CONCLUSIONS 

If very massive stars - exceeding the currently accepted 150 Mq 
limit - were to exist in the local universe, they would be: 

(i) Located in high mass (> 10^ Mq), very young (< 2 Myr) 
star clusters; 

(ii) Visually the brightest stars in their host cluster, since L oc 



M^"^ for zero age main sequence stars above 85 M^q and surface 
temperatures remain approximately constant for the first 1.5 Myr 

(iii) Possess very powerful stellar winds, as a result of the mass 
dependence of the Eddington parameter, Fe oc L/M ^ M°^ for 
such stars. 150-3 00 M© zero age main sequence stars, for which 
Fe ~ 0.4 - 0.55, would likely possess an O supergiant morphology, 
while a Wolf-Rayet appearance would be likely to develop within 
the first 1-2 Myr, albeit with significant residual surface hydrogen. 

In this study we have presented spectroscopic analyses of 
bright WN stars located within R136 in the LMC and NGC 3603 
in the Milky Way that perfectly match such anticipated charac- 
teristics. The combination of line blanketed spectroscopic tools 
and contemporary evolutionary models reveals excellent agreement 
with dynamical mass determinations for the components of Al in 
NGC 3603, with component B possessing a higher initial mass of 
~170 M0, under the assumption that it is single. Application to 
the higher temperature (see also Riihling 2008) brighter members 
of R136 suggests still higher initial masses of 165 - 320 M©. Ow- 



^ The exponent flattens further at the highest masses, such that L oc M^-^ 
for zero age main sequence stars between 300-500 Mq. 
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Table 9. Compilation of stars within R136/30 Dor, NGC 3603 and the 
Arches cluster whose initial masses exceed ^150 Mq according to evo- 
lutionary models presented here. Photometry, extinctions and bolometric 
corrections are presented in this study, with the exception of the Arches 
cluster for which we follow Espinoza et al. (2009), Kim et al. (2006) and 
Martins et al. (2008), respectively. Known or suspected binaries are marked 
with *. 



Star 


SpType 


niKs 


Ak, 


(m-M)o Mks 


BCk, 


MboI 


Minit 






mag 


mag 


mag 


mag 


mag 


mag 


M© 


R136al 


WN5h 


11.1 


0.2 


18.45 


-7.6 


-5.0 


-12.6 


320 


R136a2 


WN5h 


11.4 


0.2 


18.45 


-7.3 


-4.9 


-12.2 


240 


R136C 


WN5h 


11.3 


0.3 


18.45 


-7.4 


-4.7 


-12.1 


220* 


Mk34 


WN5h 


11.7 


0.3 


18.45 


-7.1 


-4.8: 


-11.9: 


190* 


Arches F6 


WN8-9h 


10.1 


3.1 


14.5 


-7.5 


-4.0 


-11.5 


>185 


NGC 3603 B WN6h 


7.4 


0.6 


14.4 


-7.5 


-3.9 


-11.4 


166 


R136a3 


WN5h 


11.7 


0.2 


18.45 


-6.9 


-4.8 


-11.7 


165 


Arches F9 


WN8-9h 


10.6 


3.1 


14.5 


-7.0 


-4.4 


-11.4 


>165 


Arches F4 


WN7-8h 


10.2 


3.1 


14.5 


-7.4 


-3.9 


-11.3 


>150: 


NGC 3603 


WN6h 


8.0 


0.6 


14.4 


-7.0 


-4.2 


-11.2 


148 


Ala 


















Arches F7 


WN8-9h 


10.3 


3.1 


14.5 


-7.3 


-4.0 


-11.3 


>148: 


Arches Fl 


WN8-9h 


10.35 


3.1 


14.5 


-7.25 


-4.0 


-11.25 


>145: 



ing to their dense stellar winds, they have already lost up to 20% of 
their initial mass with the first '-^l.S Myr of their main sequence 
evolution (see also de Koter et al. 1998). The R136 WN5 stars 
are moderately hydrogen depleted, to which rotating models pro- 
vide the best agreement, whereas the NGC 3603 WN6 stars pos- 
sess normal hydrogen contents, as predicted by non-rotating mod- 
els at early phases. These differences may arise from different for- 
mation mechanisms (e.g. stellar mergers). Wolff et al. (2008) found 
a deficit of slow rotators within a sample of R136 early-type stars, 
although statistically significant results await rotational velocities 
from the VLT-FLAMES Tarantula Survey (Evans et al. 2010a). 

We have assessed the potential that each of the R136 stars 
represent unresolved equal-mass binary systems. SINFONI spec- 
troscopic observations argue against close, short period binaries 
(Schnurr et al. 2009). R136c possesses an X-ray luminosity which 
is a factor of '-^ 100 times greater than that which would be expected 
from a single star and is very likely a massive binary. In contrast. 
X-ray emission from R136a exceeds that expected from single stars 
by only a factor of ^3, to which multiple wind interactions within 
the cluster will also contribute. If the R136a stars possessed simi- 
lar X-ray properties to R136c and NGC 3603C, its X-ray emission 
would be at least a factor of 15 times higher than the observed lumi- 
nosity. At most, only one of the WN sources within R136a might be 
a very long period, large separation (^300 AU), equal-mass binary 
system if as little as 30% of the shock energy contributes to the 
X-ray luminosity. Dynamical effects would harden such systems 
on <CMyr timescales. We cannot rule out shorter-period, unequal- 
mass binary systems of course, but such cases would have little 
bearing upon our mass limit inferences. 

In Table |9] we present a compilation of stars in R136 (30 Do- 
radus), NGC 3603 and the Arches cluster whose initial masses 
challenge the currently accepted upper mass limit of ^150 M©. 
Although the formation of high mass stars remains an unsolved 
problem in astrophysics (Zinnecker & Yorke 2007), there are no 
theoretical arguments in favour of such a limit at 150 M© (e.g. 
Klapp et al. 1987) - indeed Massey & Hunter (1998) argued against 
an upper stellar limit based on R136 itself. Observations of the 



Arches cluster provide the primary evidence for such a sharp mass 
cutoff (Figer 2005). However, as Table |9] illustrates, contemporary 
photometry and foreground extinction towards the cluster coupled 
with the spectroscopy results from Martins et al. (2008) suggest 4- 
5 stars initially exceed ?^150 M©, with an estimate of >185 M© 
for the most luminous star. Recall Martins et al. (2008) used iden- 
tical spectroscopic tools to those employed in the current study, 
plus near-IR spectroscopic observations which also form a central 
component of our study. On this basis the Arches cluster would no 
longer be a statistical oddity. 

Monte Carlo simulations for various upper stellar mass lim- 
its permit estimates of the revised threshold from which ^-^300 M© 
is obtained. It may be significant that both NGC 3603 and R136 
are consistent with an identical upper limit, in spite of their differ- 
ent metallicities. Oey & Clarke (2005) obtained an upper limit of 
^500 M0, primarily from a maximum stellar mass of 120 - 200 
M© inferred for stars within R136a itself at that time. 

Would there be any impact of an upper mass limit of order 
^-^300 M© on astrophysics? Population synthesis studies are widely 
applied to star-forming regions within galaxies, for which an upper 
mass limit of '-^ 120 M© is widely adopted (Leitherer et al. 1999). 
A number of properties are obtained from such studies, including 
star formation rates, enrichment of the local interstellar medium 
(ISM) through mechanical energy through winds, chemical enrich- 
ment and ionizing fluxes. A higher stellar mass limit would increase 
their global output, especially for situations in which rapid rotation 
leads to stars remaining at high stellar temperatures. 

Could the presence of very high mass stars, such as those dis- 
cussed here, be detected in spatially unresolved star clusters? The 
high luminosities of such stars inherently leads to the development 
of very powerful stellar winds at early evolutionary phases (1-2 
Myr). Therefore the presence of such stars ensures that the inte- 
grated appearance of R136a (and the core of NGC 3603) exhibits 
broad He II A4686 emission, which at such early phases would not 
be the case for a lower stellar limit. Other high mass clusters, wit- 
nessed at a sufficiently early age, would also betray the presence of 
very massive stars through the presence of He II emission lines at 
A 1640 and A4686. Hitherto, such clusters may have been mistaken 
for older clusters exhibiting broad helium emission from classical 
Wolf-Ray et stars. 

Finally, how might such exceptionally massive stars end their 
life? This question has been addressed from a theoretical perspec- 
tive (Heger et al. 2003), suggesting Neutron Star remnants fol- 
lowing Type Ib/c core-collapse supernovae close to solar metal- 
licities, with weak supernovae and Black Hole remnants at LMC 
compositions. Rapid rotation causes evolution to proceed directly 
to the classical Wolf-Rayet phase, whereas slow rotators, such as 
the NGC 3603 WN6h stars, will likely produce a rj Car-like Lumi- 
nous Blue Variable phase. Extremely metal-deficient stars exceed- 
ing ~140 M© may end their lives prior to core-collapse (Bond et 
al. 1984). They would undergo an electron-positron pair-instability 
explosion during the advanced stages that would trigger the com- 
plete disruption of the star (Heger & Woosley 2002). 

Until recently, such events have been expected from solely 
metal-free (Population III) stars, since models involving the col- 
lapse of primordial gas clouds suggest preferentially high charac- 
teristic masses as large as several hundred M© (Bromm & Lar- 
son 2004). Langer et al. (2007) have demonstrated that local pair- 
instability supernovae could be produced either by slow rotating 
moderately metal-poor (< 1/3 Z©) yellow hypergiants with thick 
hydrogen-rich envelopes - resembling SN 2006gy (Ofek et al. 
2007; Smith et al. 2007) - or rapidly rotating very metal-deficient 
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(< 10~^ Zq) Wolf-Rayet stars. Therefore it is unlikely any of 
the stars considered here are candidate pair-instability supernovae. 
Nevertheless, the potential for stars initially exceeding 140 M© 
within metal-poor galaxies suggests that such pair-instability su- 
pernovae could occur within the local universe, as has recently 
been claimed for SN 2007bi (Gal- Yam et al. 2009, see also Langer 
2009). 

Finally, close agreement between our spectroscopically de- 
rived mass-loss rates and theoretical predictions allows us to syn- 
thesise the ZAMS appearance of very massive stars. We find that 
the very highest mass progenitors will possess an emission-line ap- 
pearance at the beginning of their main- sequence evolution due 
to their proximity to the Eddington limit. As such, spectroscopic 
dwarf 02-3 stars are not anticipated with the very highest masses. 
Of course the only direct means of establishing stellar masses is 
via close binaries. For the LMC metallicity, Massey et al. (2002) 
have obtained dynamical masses of a few systems, although spec- 
troscopic and photometric searches for other high-mass, eclipsing 
binaries in 30 Doradus are in progress through the VLT-FLAMES 
Tarantula Survey (Evans et al. 2010a) and other studies (O. Schnurr 
et al. in prep.). Nevertheless, it is unlikely that any other system 
within 30 Doradus, or indeed the entire Local Group of galaxies 
will compete in mass with the brightest components of R136 dis- 
cussed here. 
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